timer on 1
use "Data\MFPS Child Health Paper Replication Data - Child Level.dta", clear


* Set Panel Variables
xtset caseid_n year
gen health_2017 = health_visits
replace health_2017 = L.health_visits if year == 2018

* Set Seed
set seed 2611 

*Create global of age covariates
mkspline age1 6 age2 12 age3 18 age4 24 age5 30 age6 = child_age_months


/*create dummy variables for missing data
local cov work_bl ever_use_bl curr_use_bl age_bl edu_primary_bl sex_age_bl birth_order religion_r_bl ethnicity_r_bl total_birth_bl child_sex child_age_months chi1m health_visits_bl wealthscore_bl
foreach i of local cov{ 
	gen `i'_missing = 0 
	replace `i'_missing = 1 if `i'==.
	replace `i' = 0 if `i'==.
}
*/
*Generate month of birth and cluster FE
tab chi1m, gen(MOB)
tab cluster_num, gen(clust)

*Generate panel variables for wave specific variables
gen next_preg_wave = 0 
replace next_preg_wave = next_preg_2017 if year ==2017
replace next_preg_wave = next_preg_2018 if year == 2018
label var next_preg_wave "Birth Spacing"

gen trans_count_wave = 0 
replace trans_count_wave = w2_trans_count if year == 2017
replace trans_count_wave = w3_trans_count if year == 2018 
label var trans_count_wave "Transportation Count"

gen couns_count_wave = 0 
replace couns_count_wave = w2_couns_count if year == 2017
replace couns_count_wave = w3_couns_count if year == 2018 
label var couns_count_wave "Counseling Count"

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
* Mediation analysis assuming independence between mediators
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

	
	////////////////////////////////////////////////////////////////////////////
	*Set globals
	** Need to manually drop colinear variables and perfect predictors because
	*** medeff will not do this automatically and the matrix will not be PSD
	////////////////////////////////////////////////////////////////////////////
		global agevars1 = "MOB1 MOB3-MOB11 age2-age4"
		global agevars2 = "MOB1 MOB3-MOB11 age4-age6"
		
		global X1 = "" //Need to manually drop colinear or missing coef
		global X2 = "clust1-clust4 child_sex ever_use_bl curr_use_bl age_bl sex_age_bl birth_order edu_primary_bl total_alive_bl religion_r_bl ethnicity_r_bl"
		
	local mediators health_2017 next_preg_wave trans_count_wave w13_child_hrs curr_use
	//////////////////////////////////////////////////
	* Set up results export 
	//////////////////////////////////////////////////

	*Format Table Header
	putexcel set "Results\Mediation\mediation effects$name.xlsx", replace
		putexcel A1 = ("Mediation Results")
		putexcel B2 = ("ACME")
		putexcel C2:D2 = ("95% Confidence Invervals")
		putexcel E2 = ("$\rho$ where ACME = 0")
		putexcel F2 = ("ITT on mediator: Beta")
		putexcel G2 = ("ITT SE")
		putexcel H2 = ("Observations")
		
	//////////////////////////////////////////////////
	* Mediate HAZ results 
	//////////////////////////////////////////////////
	local j = 3
	foreach med of local mediators{
		local name: variable label `med'
	
		forvalues i = 1/2 {

			* LSEM Mediation
			medeff (regress `med' treatment $agevars1 ${X`i'}) (regress haz06 treatment `med' $agevars1 ${X`i'}) if year==2017 & ///
			index_dummy==1 & born_during_study==0, mediate(`med') treat(treatment) vce(robust) sims(2000)
			*Format results into a table
			putexcel A`j' = ("Outcome: HAZ, Mediator: `name', Covariate Set: `i'")
				putexcel B`j' = (r(delta1))
				putexcel C`j' = (r(delta1lo))
				putexcel D`j' = (r(delta1hi))

		*Sensitivity Analysis
		medsens (regress `med' treatment $agevars1 ${X`i'}) (regress haz06 treatment `med' $agevars1 ${X`i'}) if year==2017 & ///
		index_dummy==1 & born_during_study==0, mediate(`med') treat(treatment) sims(2000)
				putexcel E`j' = (r(errcr))
			
		*ITT Effect 
		reg `med' treatment $agevars1 ${X`i'} if year==2017 & index_dummy==1 & born_during_study==0 & haz06!= ., robust
			putexcel F`j' = (_b[treatment])
			putexcel G`j' = (_se[treatment])
			putexcel H`j' = (e(N))
		
		local j = `++j'
	}

	/*Plot Sensitivity Analysis for Covariate set 3
	drop _med_updelta0 _med_lodelta0 _med_rho _med_delta0 _med_rho
	medsens (reg `med' treatment $agevars1 ${X3}) (regress haz06 treatment `med' $agevars1 ${X3}) if year==2017 & index_dummy==1 & ///
		born_during_study==0, mediate(health_visits) treat(treatment) sims(500)
	
	twoway rarea _med_updelta0 _med_lodelta0 _med_rho, bcolor(gs14) || line _med_delta0 _med_rho, lcolor(black) ytitle("ACME", size(small)) ///
	title("Outcome: HAZ; Mediator: `name'", size(small)) xtitle("Sensitivity parameter: {&rho}") legend(off) scheme(sj)

	graph save a, replace
	*/
	//////////////////////////////////////////////////
	* Mediate CREDI results with health visits in 2017
	//////////////////////////////////////////////////
	forvalues i = 1/2 {

		* LSEM Mediation
		medeff (regress `med' treatment $agevars2 ${X`i'}) (regress credi_z_score treatment `med' $agevars2 ${X`i'}) if year==2018 & ///
		index_dummy==1 & born_during_study==0, mediate(`med') treat(treatment) vce(robust) sims(2000)
		*Format results into a table
			putexcel A`j' = ("Outcome: CREDI Z-Score, Mediator: `name', Covariate Set: `i'")
				putexcel B`j' = (r(delta1))
				putexcel C`j' = (r(delta1lo))
				putexcel D`j' = (r(delta1hi))

		*Sensitivity Analysis
		medsens (regress `med' treatment $agevars2 ${X`i'}) (regress credi_z_score treatment `med' $agevars2 ${X`i'}) if year==2018 & ///
		index_dummy==1 & born_during_study==0, mediate(`med') treat(treatment) sims(2000)
			putexcel E`j' = (r(errcr))
			
		*ITT Effect 
		reg `med' treatment $agevars2 ${X`i'} if year==2018 & index_dummy==1 & born_during_study==0 &  credi_z_score!= ., robust
			putexcel F`j' = (_b[treatment])
			putexcel G`j' = (_se[treatment])
			putexcel H`j' = (e(N))
		
		local j = `++j'
	}

	/*Plot Sensitivity Analysis for Covariate set 3
	drop _med_updelta0 _med_lodelta0 _med_rho _med_delta0 _med_rho
	medsens (reg `med' treatment $agevars2 ${X3}) (regress credi_z_score treatment `med' $agevars2 ${X3}) if year==2018 & ///
	index_dummy==1 & born_during_study ==0, mediate(health_2017) treat(treatment) sims(500)
	
	twoway rarea _med_updelta0 _med_lodelta0 _med_rho, bcolor(gs14) || line _med_delta0 _med_rho, lcolor(black) ytitle("ACME", size(small)) ///
	title("Outcome: CREDI Z-Score; Mediator: Health Visits (2017)", size(small)) xtitle("Sensitivity parameter: {&rho}") legend(off) scheme(sj)

	graph save b, replace
	*/
	}
	
	

timer off 1
timer list 1